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The application of stability theory in Laminar Flow Control 
(LFC ) research requires that density and velocity profiles be 
specified throughout the viscous flow field of interest. These 
profile values must be as numerically accurate as possible and 
free of any numerically induced oscillations. To date, a high 
percentage of the boundary-layer solutions for use in three- 
dimensional (3-D) stability analysis programs have been obtained 
using quasi-3-D procedures due to the general unavailability of 
3-D boundary-layer programs. 

Several mature 3-D boundary-layer programs based on finite- 
difference techniques (note: integral approaches are not 

considered in the present paper) currently exist, for example, 
references 1 through 4. Unfortunately, these programs are not 
available for general use by the LFC research community due to 
company proprietary controls. In addition to proprietary 
controls, most of the programs are not structured or adequately 
documented for ease of use; that is, their use would generally 
require direct interaction with and/or assistance from the group 
or individual that developed the program. 

Guidelines for the present research project are presented in 
figure 1 and are as follows: (1) develop an efficient and 

accurate procedure for solving the 3-D boundary -layer equations 
for aerospace configurations; (2) develop an interface program to 
couple selected 3-D inviscid programs that span the subsonic to 
hypersonic Mach number range; and (3) document and release 
software to the LFC research community. 


• Reasons for development | 

• No general program exists that Is readily available to 
all groups involved in laminar flow control research 


• No general program exists that is coupled di recti 
the most often used inviscid programs over the su 
to hypersonic Mach number range 


y with 
Dsonlc 


• Most existing programs are at best second-order accurate 
and often degenerate to first-order accuracy 


• Guidelines for development 

• User friendly 

• High-order accuracy; G(aZ 4 >, 0<aX, ay) 2 

• Flexible interface 


• Document and release to LFC Research community 


(Finite difference numerics) 


FIGURE 1 
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ORDER OF PRESENTATION 


As a result of time constraints for the present symposium, 
it was necessary to condense the material that would have been 
included in three presentations into a simple co-authored 
paper. Consequently, no attempt will be made to either (1) give 
a review of current boundary-layer literature or (2) compare 
numerical results from the present solution procedure with 
previously published numerical and/or experimental results. The 
material to be presented is outlined in figure 2 and is divided 
into two basic categories. The interaction between these 
categories will not be discussed due to time limitations. 

Within the inviscid area the focus of the presentation will 
be on the processing of the inviscid flow field data required as 
input to the boundary-layer program. The surface Euler equations 
will be defined together with the numerical procedure for solving 
the nonlinear system. Several test cases will be presented where 
results obtained from standard interpolation procedures will be 
compared with results from the surface Euler equations. 

Within the viscous area the focus of the presentation will 
be on the numerical procedure and the establishment of its 
accuracy and numerical efficiency. The boundary-layer equations 
will not be presented in equation form; however, their 
mathematical character and the required boundary and initial 
conditions will be discussed. The reader interested in the full 
equations can find them in several publications (see refs. 1 
through 6, for example). 


• Inviscid 

• Interface program 

• Inviscid data treatment 

• Surface Euler equations 

• Boundary conditions 

• Numerical procedure 

• Test cases 

• Viscous 

• Governing equations 

• Boundary conditions 

• Initial conditions 

• Transformation 

• Solution technique 

• Numerical scheme 

• Matrix structure 

• Test cases 


FIGURE 2 
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INTERFACE PROGRAM 


The interface program provides the required linkage between 
the most appropriate inviscid code for the configuration and 
flight Mach number regime and the three-dimensional boundary- 
layer program. The functions of the interface program are 
presented in figure 3. The interface program must be capable of 
(1) generating a sufficiently refined boundary-layer grid; (2) 
calculating the metric coefficients and related geometric 
parameters for this grid; (3) interpolating the pressure field 
from the relatively coarse grid used for the inviscid flow field 
solution to the relatively fine boundary-layer grid with 
automatic smoothing capability as required; (4) solving the 
surface Euler equations on the wall boundary to obtain the 
inviscid velocity components and their gradients in the £, yi 
coordinate system; and (5) generating all required input files 
for the boundary-layer program. In addition to these primary 
functions, the interface program prepares diagnostic plot files 
and will act as the interactive link for viscous-inviscid 
interaction calculations. 

Success (measure of acceptance and use of software by LFC 
research groups) of the boundary-layer program will depend, to a 
large extent, on the structure of the interface software. 


• General function 

• Provide the required linkage between selected inviscid 
programs that span the subsonic to hypersonic Mach 
number range and the boundary- layer program under 
development 

• Specific functions 

• Generation of boundary- layer grid and related metric 
coefficients 

• Interpolation of inviscid pressure distribution. from 
the relatively coarse inviscid grid to the sufficiently 
fine boundary- layer grid 

• Solution of the surface Euler equations to obtain 
the inviscid velocity distribution over the 
boundary- layer grid 

• Generation of input and output files for the boundary- 
layer program 

• Interface for viscous-inviscid interaction 


figure 3 



TREATMENT OF INVISCID DATA 


The sensitivity of stability theory to oscillations in the 
viscous flow field requires that the inviscid data and body 
geometry (metric coefficients, etc.) be carefully treated in 
order to avoid the introduction of nonphysical oscillations. 
Alternate approaches to inviscid flow-field data treatment are 
presented in figure 4. A conventional approach has been to use 
spline interpolation to transfer data from the relatively coarse 
inviscid grid to the sufficiently fine boundary-layer grid. 

These interpolated results are then numerically smoothed. The 
interpolation approach is the most straightforward; however, 
large errors can occur in the interpolation and arbitrary 
smoothing of vector quantities. This approach yields results 
that are not consistent; that is, they do not satisfy the 
governing inviscid equations within an acceptable error bound. 

The present approach is to interpolate the pressure field 

from the coarse inviscid grid to the fine boundary-layer grid. 

This pressure distribution P ,,(£,n), together with known initial 

w 

and boundary values, is used in the numerical solution of the 
surface Euler equations. The advantage of the approach is that 
exact values of u and v can be enforced on certain boundaries, 
the velocity components (u,v) and their gradients in 
the £, n plane are smooth, and the numerical results are 
consistent . 


• conventional approach 

• Interpolate 0, v, w, P from 
inviscid solution grid to 
boundary- layer grid 


• Present approach 

• Interpolate pressure field 
from inviscid solution grid 
to boundary- layer grid 


Transform u, v, w to u, v 
on boundary-layer grid 


Surface element 


Z Y 


• Disadvantages 

• Errors from interpolation 
of vector quantities can 
be very large 



Numerically solve the surface 
Euler system to generate 
u, v on boundary- layer grid 


Inviscid grid 

Viscous grid 


Advantages 


• "Exact" values of u, v can 
be enforced on known 
boundaries 


• Difficult to obtain 
smooth distributions 

• Numerical results are 
not consistent 


• Smoothness of solution of 
u. v and their gradients 

• Numerical results are 
consistent 


FIGURE 4 


SURFACE EULER EQUATIONS 


The surface Euler equations are obtained from the three- 
dimensional boundary-layer equations in the limit as 5 -*• <*>. The 
system obtained from the £ and n momentum equations are presented 
in figure 5. The system consists of two nonlinear hyperbolic 
equations for the E, and n momentum equations and an algebraic 
relationship for the energy equation (perfect gas). The b^ 
values are known functions of the grid system. The pressure is 
known from the interpolation procedure. The unknowns are 
p, u, and v, where p, u, and v are the density, ^-velocity 
component, and n-velocity components, respectively. The 
equations are a first order, nonlinear hyperbolic system where 
the local streamlines are characteristic. The system can be 
solved with an explicit march in £ while remaining fully implicit 
in n provided the initial conditions and boundary conditions can 
be specified as indicated in figure 5. 


A method for calculating the inviscid velocity vectors from 
a specified pressure field was first reported by Cebeci and Meier 
(ref. 6) for incompressible flow around an ellipsoid at angle of 
attack. Vollmers (ref. 7) and Gleyzes and Cousteix (ref. 8) 
improved the approach by integrating the Euler equations. 

Anderson (ref. 3) described a method using experimental pressure 
distributions on rotating turbine blades using the Euler 
equations . The current procedure extends these concepts to 
complex geometry associated with wings and fuselages using a 
second-order numerical scheme. 


j 

i 



• Equation system obtained from boundary- layer equations 
In limit as — 

• Equations 

5 -momentum 


>io n u-i^ 


+ °12 v §^ + (t, 12 u ' + b 14 uv + D 15 v ' n * ,b 16it + b 17fn * 


3 P 


3 P , 


n -momentum 


3 P 


p[t) 2 iU§! + b 22 v §7i + lt) 23 u2 + D 24 uv + °25 v2)1 * ' u 265c 
i Character of system 

•First order, nonlinear hyperbolic system 
•Local streamlines are characteristics 
•Explicit march in e; implicit in n 


3 P , 
“2737 


n 
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BOUNDARY CONDITIONS FOR FUSELAGE "TYPE" GEOMETRY 


As previously noted in figure 5, initial and boundary 
conditions must be specified for physical flows corresponding to 
aerospace configurations. These conditions are presented in 
figure 6(a) for a fuselage type body. The fuselage type body 
with a plane of symmetry relative to the flow field has well 
defined initial and boundary conditions. The initial conditions 
at the stagnation point in the physical plane are u = v = 0. 

These conditions from the required initial values along 

the ^ line in transformed space. Along both the leeward and 

windward lines of the symmetry plane in physical space one can 

specify the conditions v = 0 and 3u/3n = 0. These form the 

boundary conditions in the transformed plane for the 

lines n = n and n = ri 

o max 


• Physical space (X,Y,Z) • Transformed space (c, n) 




v = 0; frr = 0 


FIGURE 6 (A) 
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BOUNDARY CONDITIONS FOR WING GEOMETRY 


Boundary and initial conditions for wing type geometry are 
somewhat more complex than those for fuselages. A schematic of a 
typical wing element is presented in figure 6(b) where a and A 
denote the angle of attack and leading-edge sweep, 
respectively. At angle of attack, the leading-edge attachment 
line will be displaced from the geometrical leading edge as shown 
in the sketch. The location of the attachment line is not known 
a priori and must be obtained by an iterative solution procedure 
that positions the line such that u = 0 to within a specified 
error bound. This becomes the % = 5 line in the transformed 

plane. The remaining initial condition along the £ = 5 line is 

obtained from the inviscid flow field results (see ref. 9). The 
boundary conditions along either the root chord or symmetry plane 
chord and a chord near, but not located at the wing tip, must be 
specified. The extrapolation conditions shown in the 5 , n plane 
in figure 6(b) for the wing tip and root chord regions have been 
used for the results contained in the present material. A more 
detailed discussion of the initial conditions for the leading- 
edge attachment line and boundary condition studies for the wing 
tip and root regions is presented in reference 9. 


• Physical space (X,Y,Z) • Transformed space ( e,n j 



FIGURE 6 (B) 



NUMERICAL PROCEDURE 


As previously mentioned, the surface Euler equations are 
hyperbolic. In order to satisfy the stability requirement, the 
difference stencil must be rotated in relation to the cross flow 
velocity as shown in figure 7(a). 


Difference stencil 
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FIGURE 7 (A) 
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DIFFERENCE APPROXIMATIONS 


The point upwind differences are used to approximate 
derivatives with respect to Z (see fig. 7(b)). For derivatives 
with respect to n a 5-point stencil is used and rotated to 
satisfy stability where e = 1 for v < 0; e = -1 for v > 0. 

Nonlinear terms are treated in an iterative sense by lagging 
the nonlinear quantity by one iterative step. The final set of 
discretized equations are "locally block diagonal." 
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FIGURE 7(B) 
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PRESSURE COEFFICIENT CONTOURS 


O RIGINAL -TAGS- 1 5 

OE -POOR QUALIFY 


A test case of interest to the Langley Research Center's LFC 
program involves experimental and analytical studies of the 
stability characteristics of a fuselage forebody shown in figure 
8. Figure 8(a) presents contours of the pressure coefficient 
distribution interpolated from a relatively coarse inviscid 
solution onto a fine boundary-layer grid. The initial inviscid 
distribution was obtained using the method of reference 10. 






INVISCID FLOW-FIELD CONTOURS 


The pressure contours in figure 8(a) indicate a smooth and 
"well behaved" flow field; however, close examination of the 
cross flow velocity contours in figure 8(b) indicates a physical 
problem not easily solved by traditional boundary-layer 
procedures. The solid line on figure 8(b) indicates the locus of 
a v = 0 contour line. In the region above this line v > 0 
(positive cross flow) while below the line v < (negative cross 
flow) . Consequently, both the leeward and windward lines of the 
plane of symmetry are inflow lines. This presents a severe 
problem for traditional boundary-layer approaches where one 
normally obtains a numerical solution of the reduced set of 
boundary-layer equations along the windward line (generally an 
outflow line). Using this solution as one of the two required 
initial data planes, the solution can be marched in the direction 
of increasing n (direction of positive cross flow velocity) to 
obtain the solution for £ = n D < n < n m ax' an< -* ^ * T ^ e 

surface Euler equations could be numerically integrated 

for E=E.,n < n < n because the entire n~line was treated 
s o max 

implicitly. 



FIGURE 8(B) 
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SURFACE EULER VERSUS STANDARD INTERPOLATION 


Comparisons of interpolated u and v values with those 
obtained from solving the surface Euler equations are presented 
in figure 9. The interpolated values were carefully smoothed to 
get the results presented in figure 9(a) as opposed to no 
smoothing for the velocity field obtained from the surface Euler 
equations. In addition, the interpolated values are not 
consistent . 


Cessna fuselage forebody 

= 0,3; a = o° 




FIGURE 9 (A) 
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CROSS FLOW VELOCITY CONTOURS 


The carpet plots of cross flow velocity presented in figure 
9(b) are indicative of oscillations resulting from interpolation 
of the velocity components. Stability theory is critically 
sensitive to oscillations of cross flow velocity? consequently, 
the surface Euler approach should be used to obtain smooth and 
consistent input data for the boundary-layer solution. 


Cessna fuselage forebody 

M°° - 0 . 3 j cx = 0 ° 



FIGURE 9(B) 


COMPARISON OF SURFACE EULER WITH INTERPOLATION 


Numerical results are presented in figure 10 for a wing 
having an NACA-0012 airfoil section. The inviscid flow field was 
obtained using the method of reference 10. Comparisons of the 
interpolated and smoothed streamwise velocity are compared with 
surface Euler results in figure 10(a). The agreement in 
magnitude appears to be satisfactory with the exception of a 
region just downstream of the leading-edge attachment line. The 
agreement appears to be independent of span location. 


NACA-0012) yv = 30°; a. = 8°; = 0.22; AR = 2.1 


(upper surface) 



X 


FIGURE 10(A) 
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COMPARISON OP SURFACE EULER WITH INTERPOLATION 


The stability equations will be integrated forward in space 
beginning at the leading-edge attachment line. A closer 
examination of the results presented in figure 10(a) is presented 
as a derivative with respect to x in figure 10(b). The 
oscillation in 3U/35 presented in figure 10(b) for the 
interpolated values is not acceptable. It can be seen that the 
surface Euler eguations yield smooth values of 30/3£ without 
reguiring smoothing. These values are consistent with the 
governing eguations. 1 

I 


NACA-0012; = 30°j c* = 8°j = 0.22; AR = 2.1 

(upper surface) 

A 



0 .025 .050 .075 


FIGURE 10(B) e 
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THREE-DIMENSIONAL BOUNDARY-LAYER SYSTEM 


Researchers involved in the numerical aspects of LFC are 
familiar with the three-dimensional boundary-layer equations; 
consequently, the equations are not presented. Readers 
interested in the equations are refered to references 1 through 
5. Certain characteristics of the system are presented in figure 
11 and warrant further discussion. 

The system consists of four nonlinear partial-differential 
equations (PDEs) of mixed type. The system is hyperbolic in 
planes parallel to the boundary surface (£,n planes) and 
parabolic in the direction normal to these planes. 

For perfect gas flow, the system is closed with algebraic 
equations for state and viscosity. The unknowns are u, v, w, T 
and p representing the velocity components in the 
directions, static temperature and density, respectively. 

Pressure is constant normal to the wall boundary. 

The present work focuses on laminar flow; however, the 
software has been structured to model the turbulent terms using 
eddy viscosity and eddy conductivity closure. Careful studies 
will be made in this area at some future date. This would allow 
a uniform approach for boundary-layer calculations to be 
developed for laminar, transitional, and turbulent flows. 


• Equations 

• Continuity 

• Momentum (2; X, Y) 

• Energy 

• State; viscosity } 


Nonlinear PDE 
Algebraic 


• Turbulence 

•Eddy viscosity/eddy conductivity j Not LFC requirement 


FIGURE 11(A) 


THREE-DIMENSIONAL BOUNDARY-LAYER SYSTEM 


Boundary conditions have been treated as generally as 
possible. Wall mass transfer can be specified together with 
either wall temperature or heat transfer distributions. Edge 
values are obtained from the interface software solutions of the 
surface Euler equations. Initial values can currently be 
generated for fuselage and wing geometries. The reader 
interested in the special equation sets for the initial value 
planes is referred to reference 5. 


• Boundary conditions 
• Wall boundary 
• U; V (no slip) 

•W w = W w (X,Y) 

•T w = T w (X,Y) 
or q w = q w (X,Y) 


• Edge values 

•Surface Euler system 


# Initial values 

• Fuselage "type" geometry • Wing "type" geometry 

• 3-D stagnation point •Leading edge attachment line 

• Symmetry plane •Symmetry plane 

•Root chord approximation 

FIGURE 11(B) 
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NUMERICAL SCHEME 


The difference stencil together with a sketch depicting the 
Raetz influence principle (ref. 5), and the difference 
approximations are presented in figure 12(a). Two separate 
stencils are used for convective derivatives in the £,n-plane. 
The particular stencil utilized at each £ point depends upon the 
sign of the cross flow velocity. For negative cross flow the 
Krause et al. (ref. 11) zig-zag scheme is used to satisfy the 
zone of influence principle. 


• Difference stencil 
Positive cross flow 



Negative cross flow 

1-1J+1 


LLk-1/2 
O i,j,k-l 



k-1/2 planeljjL 
ij-l 


• Raetz influence principle 

Domain of dependence -* — Domain of influence 

Initial data 



Characteristic surfaces 


FIGURE 12 (A) 
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DIFFERENCE APPROXIMATIONS 


The difference approximations are presented in figure 
12(b). The approximations are second-order accurate in planes 
parallel to the wall and fourth-order accurate in the direction 
normal to the wall boundary. In the planes parallel to the wall, 
a weighted three-point upwind scheme is used to approximate the 
convective partial-derivatives (see ref. 12). 

In the direction normal to the wall boundary, a 2-point 
compact scheme with fourth-order accuracy is used (see refs. 13 
and 14). 


>e - derivatives 
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NUMERICAL SCHEME 

Difference relations 
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FIGURE 12(B) 



SOLUTION TECHNIQUE 


Substitution of the difference relations into the 
transformed governing eguations yields a set of seven nonlinear, 
coupled difference eguations in seven unknowns (see fig. 13). 
Newton's method is used to linearize the system. The linearized 
system is solved using an efficient block-tridiagonal matrix 
inversion based on LU factorization. 

The two momentum eguations yield a 4x4 block tridiagonal 
system in F, F', G, G', where F = u/ue , G = v/ve , F' = 3F/3e and 
G’ = 3G/3 The energy eguation yields a 2x2 block for H and 
H'. The normal component of velocity, w is updated from the 
finite-difference form of the continuity eguation. The system is 
solved in an iterative loop until convergence is obtained within 
a fixed error bound. 


• Difference equations are coupled and nonlinear 

• Newton's method used for linearization 

Q p = Q P_1 + <sQ P_1 

• The linearized system of equations in <sq forms a block 
tri-diagonal system of the form 

lAlsQ^ 1 + |B(sQ k + )Cl«Q k+1 = R k 

Where 

lA( (c( are block matrices of system 

Q is solution vector 

4x4 block tri-diagonal for momentum eqs 

2x2 block tri-diagonal for energy eq 

w updated from difference form of continuity eq 


FIGURE 13 
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ACCURACY OF METHOD 


The results of a study of the accuracy of the numerical 

scheme together with the effect of the number of grid points 

normal to the wall boundary for a fixed c are presented in 

max 

figure 14. The truncation error study for Blasius flow clearly 
indicates a slope of 4 to 1 as compared with 2 to 1 for the 
second -order scheme of reference 12. 

The numerical advantage of the fourth-order method is shown ( 

for 3-D stagnation point flow where streamwise velocity profiles ! 

( u / Ug ) are presented for a fixed S max for 40, 16, and 8 points j 

normal to the wall boundary (note that only a sample of the > 

computed points is presented for clarity). Comparisons of the 
results for NPZ = 8 with NPZ = 40 indicate that the correct 
profile values can be predicted with as few as 8 points. This 
comparison represents a speed/storage advantage of 5 to 1. 


Log 10' 

(error 

Of Cj) 


FIGURE 14 


Blasius flow 3-D stagnation point flow 


■1-7 r — Fourth -order scheme 

I ^ _ — — Cornnrl -nrrlor crhomo 



Log 10 , (no. of points) U/U E 
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CYLINDER ON PLAT PLATE 


Numerical results for incompressible flow past a cylinder 
mounted normal to a flat plate are presented in figure 15. Total 
and cross flow velocity profiles are presented for 50, 25, and 12 
points normal to the wall as well as compared with the second- 
order theory of reference 12. The present fourth-order method 
produces excellent results for as few as 12 grid points normal to 
the wall boundary for this test problem. 


Total velocity 


Cross flow velocity 


Zeta 



FIGURE 15 
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PROLATE SPHEROID AT ANGLE OF ATTACK 


Numerical results for F* and G' w for incompressible flow 

past a prolate spheroid at 8° angle of attack is presented in 
figure 16. The results are presented for three x locations where 
n varies from 0° on the most windward line of symmetry to 180° on 
the leeward line. Of particular importance is the reverse region 
of cross flow: G' w < 0; x = 1.52; 60° < n < 180° (note that 

G’ w = (3G/3 C ) w ; G = v/v 0 ). 


F i 


Z 



FIGURE 16 



SWEPT WING FLOW 


The final swept wing flow test case presented is for a wing 
having an NACA-0012 airfoil section. A schematic of the geometry 
and test conditions is presented in figure 17(a) together with 
calculated values of the chordwise ( 5 ) and spanwise ( n ) skin 
friction coefficients. 


Schematic of geometry 

= 0.22 

N r = 7.7 - 10 6 

cx = 8° 

a = 30 ° 



Skin friction coefficient 
4 - 10' 8 



0 .2 .4 .6 .8 1.0 


FIGURE 17(A) 
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SWEPT WING FLOW: VELOCITY PROFILES 


Calculated total and cross flow velocity profiles are 
presented in figure 17(b). The numerical results were obtained 
for grid point distributions normal to the wall boundary of 12, 
20, and 50 points. The profile calculations are presented at 
chord stations of 0, 50, and 90 percent. The present fourth- 

order method obtained accurate results with as few as 12 points 
normal to the wall boundary. 


Total velocity 

X/C = .92 40)0.0 

35 U 


»r! 0 presem 
scheme 
— NPZ = 50 
1-0 20 

□ 12 


Cross flow velocity 
-4 


z/c 


X/C = .92 


= .50 



,2 .4 .6 .8 1.0 1.2 

UT/UTE 


-3 -2 



FIGURE 17(B) 



CONCLUDING REMARKS 


In concluding, the current software package (interface 
software; boundary-layer software) is operational and has been 
tested for several 3-D flows (see figure 18). 

The interface program has been found to be a dependable 
approach for developing a user friendly procedure for generating 
the boundary-layer grid and transforming an inviscid solution 
from a relatively coarse grid to a sufficiently fine boundary- 
layer grid. The surface Euler equations used for this procedure 
yield smooth inviscid velocity components and gradients along the 
boundary-layer coordinate lines. These data are consistent with 
the governing equations in satisfying the boundary-layer 
equations in the limit as the distance normal to the wall 
boundary becomes very large. The interface program will 
eventually function as the iterative link between the selected 
inviscid software and the boundary-layer software for 
inviscid/viscous interactions studies . 

The boundary —layer program has been shown to be fourth— order 
accurate in the direction normal to the wall boundary and second- 
order accurate in planes parallel to the boundary. The fourth- 
order accuracy allows accurate calculations with as few as one- 
fifth the number of grid points required for conventional second- 
order schemes. 


• Basic programs have been developed and verified for 
several test cases 

•Surface Euler equation approach yields smooth 
and consistent inviscid edge data 

• Boundary-layer procedure has been verified as 
computationally efficient with fourth order 
accuracy 


FIGURE 18 
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